﻿using System;
using System.Linq;
using System.Collections.Generic;
using System.Collections;
using Bling.Util;
using Bling.DSL;
using Bling.Core;
using Bling.Vecs;

namespace Bling.Core {
  using Bling.Matrices;
  public partial class FlexDoubleBl<D> : DoubleBl<D, FlexDoubleBl<D>> where D : Dim<D> {
  }


  public abstract partial class FixDoubleBl<D, BRAND> : DoubleBl<D, BRAND>
    where D : Dim1D<D>
    where BRAND : FixDoubleBl<D, BRAND> {
    public static explicit operator MatrixBl<D1, D>(FixDoubleBl<D, BRAND> b) { return VecMatrix<D>.FromVecA.Instance.Make(b); }
    public static implicit operator MatrixBl<D, D1>(FixDoubleBl<D, BRAND> b) { return VecMatrix<D>.FromVecB.Instance.Make(b); }
    public static BRAND Convert(MatrixBl<D1, D> b) { 
      return ToBrandD(VecMatrix<D>.FromMatrixA.Instance.Make(b)); 
    }
    public static BRAND Convert(MatrixBl<D,D1> b) {
      return ToBrandD(VecMatrix<D>.FromMatrixB.Instance.Make(b));
    }
  }
  public partial class DoubleBl : BaseDoubleBl {
    public static implicit operator DoubleBl(MatrixBl<D1, D1> M) { return M[0, 0]; }
    public static implicit operator MatrixBl<D1, D1>(DoubleBl d) { return new MatrixBl<D1,D1>(d); }
  }

  public partial class Double2Bl : PointDoubleBl<D2, Double2Bl>, IHasXY<double, DoubleBl> {
    public static explicit operator Double2Bl(MatrixBl<D1, D2> D) { return VecMatrix<D2>.FromMatrixA.Instance.Make(D); }
    public static implicit operator Double2Bl(MatrixBl<D2, D1> D) { return VecMatrix<D2>.FromMatrixB.Instance.Make(D); }
    public MatrixBl<D1, D2> Transpose { get { return (MatrixBl<D1, D2>)this; } }
  }
  public partial class Double3Bl : PointDoubleBl<D3, Double3Bl>, IHasXYZ<double, DoubleBl> {
    public static explicit operator Double3Bl(MatrixBl<D1, D3> D) { return VecMatrix<D3>.FromMatrixA.Instance.Make(D); }
    public static implicit operator Double3Bl(MatrixBl<D3, D1> D) { return VecMatrix<D3>.FromMatrixB.Instance.Make(D); }
    public static explicit operator Double3Bl(Double2Bl P) { return new Double3Bl(P, 0d); }
    /*
    public partial class To4Op : UnaryOperatorS<IVec<double, D4>, To4Op>, INoLinqOperator, INoShaderOperator {
      public new static readonly To4Op Instance = new To4Op();
      private To4Op() { }
      public override string Name {
        get { return "to4"; }
      }
      public override Expr<IVec<double, D4>> Expand(Expr<IVec<double, D3>> argA) {
        return new Double4Bl(ToBrandD(argA), 1d);
      }
      public override Expr<IVec<double, D3>> Inverse(Expr<IVec<double, D3>> argA, Expr<IVec<double, D4>> result) {
        return (Double3Bl)((Double4Bl)result);
      }
    }
     */
    public static implicit operator Double4Bl(Double3Bl Value) {
      return new Double4Bl(Value, 1d);
      //return To4Op.Instance.Make(Value); 
    }
    public MatrixBl<D1, D3> Transpose { get { return (MatrixBl<D1, D3>)this; } }
  }
  public partial class Double4Bl : FixDoubleBl<D4, Double4Bl>, IHasXYZW<double, DoubleBl> {
    public static explicit operator Double4Bl(MatrixBl<D1, D4> D) { return VecMatrix<D4>.FromMatrixA.Instance.Make(D); }
    public static implicit operator Double4Bl(MatrixBl<D4, D1> D) { return VecMatrix<D4>.FromMatrixB.Instance.Make(D); }

    /*
    public partial class To3Op : UnaryOperatorS<IVec<double, D3>, To3Op>, INoLinqOperator, INoShaderOperator {
      public new static readonly To3Op Instance = new To3Op();
      private To3Op() { }
      public override string Name {
        get { return "to3"; }
      }
      public override Expr<IVec<double, D3>> Expand(Expr<IVec<double, D4>> argA) {
        return ToBrandD(argA).XYZ() / ToBrandD(argA).W;
      }
      public override Expr<IVec<double, D4>> Inverse(Expr<IVec<double, D4>> argA, Expr<IVec<double, D3>> result) {
        return (Double4Bl) Double3Bl.ToBrandD(result);
      }
    }*/
    public static explicit operator Double3Bl(Double4Bl Value) {
      return Value.XYZ() / Value.W;  
      //return (To3Op.Instance.Make(Value)); 
    }
    public MatrixBl<D1, D4> Transpose { get { return (MatrixBl<D1, D4>)this; } }
  }

}

namespace Bling.Matrices {
  using Bling.Core;
  public interface IMatrix<BRAND> where BRAND : IMatrix<BRAND> { }
  public interface IMatrix<D, BRAND> : IMatrix<BRAND> where BRAND : IMatrix<D, BRAND> where D : Dim2D<D> { }
  public abstract partial class BaseMatrixBl<D, BRAND> : Core.ValueBrand<IVec<double,D>, BRAND, double, D, DoubleBl>, IMatrix<D,BRAND>
    where D : Dim2D<D>
    where BRAND : BaseMatrixBl<D, BRAND> {
    public BaseMatrixBl(Expr<IVec<double, D>> Underlying) : base(Underlying) { }
    public BaseMatrixBl(DoubleBl[] Es) : base(Es) { }
    public BaseMatrixBl(DoubleBl[,] Es) : base(Convert(Es)) { }
    public BaseMatrixBl() : base() { }
    private static DoubleBl[] Convert(DoubleBl[,] Es) {
      var Es0 = new DoubleBl[Dim<D>.Count];
      for (int m = 0; m < Es.GetLength(0); m++)
        for (int n = 0; n < Es.GetLength(1); n++) {
          Es0[m * Es.GetLength(1) + n] = Es[m, n];
        }
      return Es0;
    }
    public static explicit operator FlexDoubleBl<D>(BaseMatrixBl<D, BRAND> M) { return new FlexDoubleBl<D>(M.UnderlyingD); }
  }
  public interface IMatrix<ROW,COL, BRAND> : IMatrix<Dim<ROW,COL>,BRAND>
    where BRAND : IMatrix<ROW,COL, BRAND>
    where ROW : Dim1D<ROW> where COL : Dim1D<COL> { }
  public abstract partial class BaseMatrixBl<ROW, COL, BRAND> : BaseMatrixBl<Dim<ROW, COL>, BRAND>, IMatrix<ROW,COL,BRAND>
    where ROW : Dim1D<ROW> where COL : Dim1D<COL>
    where BRAND : BaseBaseMatrixBl<ROW, COL, BRAND> {
    public BaseMatrixBl(Expr<IVec<double, Dim<ROW,COL>>> Underlying) : base(Underlying) { }
    public BaseMatrixBl(DoubleBl[] Es) : base(Es) { }
    public BaseMatrixBl(DoubleBl[,] Es) : base(Es) { }
    public BaseMatrixBl() : base() { }
    public BaseMatrixBl(Func<int, int, DoubleBl> F) : this(Convert(F)) {
      if (typeof(COL) == typeof(D1)) {
        true.Assert();
      } else if (typeof(ROW) == typeof(D1)) {
        true.Assert();
      }
    
    }
    public BaseMatrixBl(params MatrixBl<D1, COL>[] Rows) : this(Convert(Rows)) { }
    public BaseMatrixBl(params MatrixBl<ROW, D1>[] Columns) : this(Convert(Columns)) { }

    public static readonly Func<Func<int, int, DoubleBl>, BRAND> ToBrandF;

    static BaseMatrixBl() {
      var mthd = typeof(BRAND).GetConstructor(new Type[] { typeof(Func<int, int, DoubleBl>) });
      (mthd != null).Assert();
      ToBrandF = f => (BRAND) mthd.Invoke(new object[] { f });
    }
    public static implicit operator MatrixBl<ROW, COL>(BaseMatrixBl<ROW, COL, BRAND> m) { return m.Provider; }
    public static BRAND operator /(BaseMatrixBl<ROW, COL, BRAND> Matrix, DoubleBl d) {
      var v = new FlexDoubleBl<Dim<ROW, COL>>(Matrix.Provider);
      v = v / FlexDoubleBl<Dim<ROW, COL>>.Repeat(d);
      return ToBrandD(v.Provider);
    }
    public static BRAND operator *(BaseMatrixBl<ROW, COL, BRAND> Matrix, DoubleBl d) {
      var v = new FlexDoubleBl<Dim<ROW, COL>>(Matrix.Provider);
      v = v * FlexDoubleBl<Dim<ROW, COL>>.Repeat(d);
      return ToBrandD(v.Provider);
    }
    public Core.StringBl ToStringBl(DoubleBl D) {
      var v = new FlexDoubleBl<Dim<ROW, COL>>(Provider);
      return v.ToStringBl(D);
    }
    public DoubleBl this[IntBl row, IntBl col] { get { return Access2DExpr<double, Dim<ROW, COL>>.Make(Provider, row, col); } }
    public class MSet {
      internal BaseMatrixBl<ROW, COL, BRAND> Outer;
      /// <summary>
      /// Specify base-1 index.
      /// </summary>
      public DoubleBl this[IntBl row, IntBl col] {
        get { return Outer[row - 1, col - 1]; }
      }
    }
    /// <summary>
    /// For one-based based access
    /// </summary>
    public MSet M { get { return new MSet() { Outer = this }; } }


    public partial class GetRowOp : BinaryOperator<IVec<double, Dim<ROW, COL>>, int, IVec<double, COL>, GetRowOp>, INoLinqOperator {
      public new static readonly GetRowOp Instance = new GetRowOp();
      private GetRowOp() { }
      public override string Name { get { return "row"; } }
      public override Expr<IVec<double, COL>> Expand(Expr<IVec<double, Dim<ROW, COL>>> argA, Expr<int> argB) {
        var Es = new Expr<double>[Dim<COL>.Count];
        for (int i = 0; i < Es.Length; i++) Es[i] = ToBrand(argA)[argB, i];
        return Composite<double, COL>.Make(Es);
      }
    }
    public partial class GetColOp : BinaryOperator<IVec<double, Dim<ROW, COL>>, int, IVec<double, ROW>, GetColOp>, INoLinqOperator {
      public new static readonly GetColOp Instance = new GetColOp();
      private GetColOp() { }
      public override string Name { get { return "col"; } }
      public override Expr<IVec<double, ROW>> Expand(Expr<IVec<double, Dim<ROW, COL>>> argA, Expr<int> argB) {
        var Es = new Expr<double>[Dim<ROW>.Count];
        for (int i = 0; i < Es.Length; i++) Es[i] = ToBrand(argA)[i, argB];
        return Composite<double, ROW>.Make(Es);
      }
    }
    public class RowColSet<ROW0,COL0> : IEnumerable<MatrixBl<ROW0,COL0>> where ROW0 : Dim1D<ROW0> where COL0 : Dim1D<COL0> {
      internal BRAND Outer;
      internal Func<IntBl, MatrixBl<ROW0, COL0>> AccessF;
      /// <summary>
      /// Access row/column by index.
      /// </summary>
      public MatrixBl<ROW0, COL0> this[IntBl Index] { get { return AccessF(Index); } }
      public IEnumerator<MatrixBl<ROW0, COL0>> GetEnumerator() {
        var ret = new MatrixBl<ROW0, COL0>[Math.Max(Dim<ROW0>.Count, Dim<COL0>.Count)];
        for (int i = 0; i < ret.Length; i++) ret[i] = this[i];
        IEnumerable<MatrixBl<ROW0, COL0>> ret0 = ret;
        return ret0.GetEnumerator();
      }
      System.Collections.IEnumerator System.Collections.IEnumerable.GetEnumerator() { return GetEnumerator(); }
    }


    /// <summary>
    /// Access matrix rows.
    /// </summary>
    public RowColSet<D1, COL> Rows {
      get {
        return new RowColSet<D1, COL>() {
          Outer = this,
          AccessF = (row) => {
            var d = GetRowOp.Instance.Make(this, row);
            return VecMatrix<COL>.FromVecA.Instance.Make(d);
          },
        };
      }
    }
    /// <summary>
    /// Access matrix column.
    /// </summary>
    public RowColSet<ROW, D1> Columns {
      get {
        return new RowColSet<ROW, D1>() {
          Outer = this,
          AccessF = (col) => {
            var d = GetColOp.Instance.Make(this, col);
            return VecMatrix<ROW>.FromVecB.Instance.Make(d);
          },
        };
      }
    }

    private static DoubleBl[] Convert(MatrixBl<D1, COL>[] Rows) {
      (Rows.Length == Dim<ROW>.Count).Assert();
      var Ds = new DoubleBl[Dim<COL>.Count * Dim<ROW>.Count];
      for (int i = 0; i < Dim<ROW>.Count; i++) {
        for (int j = 0; j < Dim<COL>.Count; j++) {
          Ds[i * Dim<COL>.Count + j] = Rows[i][j];
        }
      }
      return Ds;
    }
    private static DoubleBl[] Convert(MatrixBl<ROW,D1>[] Columns) {
      (Columns.Length == Dim<COL>.Count).Assert();
      var Ds = new DoubleBl[Dim<COL>.Count * Dim<ROW>.Count];
      for (int i = 0; i < Dim<ROW>.Count; i++) {
        for (int j = 0; j < Dim<COL>.Count; j++) {
          Ds[i * Dim<COL>.Count + j] = Columns[j][i];
        }
      }
      return Ds;
    }

    private static DoubleBl[] Convert(Func<int, int, DoubleBl> F) {
      var Es0 = new DoubleBl[Dim<ROW,COL>.Count];
      for (int m = 0; m < Dim1D<ROW>.Count; m++)
        for (int n = 0; n < Dim1D<COL>.Count; n++) {
          Es0[m * Dim1D<COL>.Count + n] = F(m,n);
        }
      return Es0;
    }

    protected class TransposeOp : UnaryOperatorS<IVec<double, Dim<COL, ROW>>, TransposeOp>, //Ops.IMemoOperator<IVec<double, Dim<COL, ROW>>>, 
      ICallOperator {
      public new static TransposeOp Instance = new TransposeOp();
      private TransposeOp() { }
      public override string Name { get { return "Transpose"; } }
      public override Expr<IVec<double, Dim<COL, ROW>>> Expand(Expr<IVec<double, Dim<ROW, COL>>> argA) {
        return new MatrixBl<COL, ROW>((m, n) => ToBrandD(argA)[n, m]).Provider;
      }
      public override Expr<IVec<double, Dim<ROW, COL>>> Inverse(Expr<IVec<double, Dim<ROW, COL>>> argA, Expr<IVec<double, Dim<COL, ROW>>> result) {
        return MatrixBl<COL,ROW>.ToBrandD(result).Transpose;
      }
      public override Func<Microsoft.Linq.Expressions.Expression, Microsoft.Linq.Expressions.Expression> Linq(Bling.Ops.Operation<IVec<double, Dim<ROW, COL>>, IVec<double, Dim<COL, ROW>>> op) {
        if (Dim<ROW>.Count == 1 || Dim<COL>.Count == 1) return e => e;
        return base.Linq(op);
      }
    }
    public virtual MatrixBl<ROW, MID> Product<MID>(MatrixBl<COL, MID> M) where MID : Dim1D<MID> {
      return ProductOp<MID>.Instance.Make(this, M);
    }
    protected class ProductOp<MID> : 
      Ops.Operator<IVec<double, Dim<ROW, COL>>, IVec<double, Dim<COL, MID>>, IVec<double, Dim<ROW, MID>>, ProductOp<MID>>,
      Ops.IMemoOperator<IVec<double, Dim<ROW, MID>>> where MID : Dim1D<MID> {

      public new static ProductOp<MID> Instance = new ProductOp<MID>();
      private ProductOp() { }

      public override Func<SmartString, SmartString, SmartString> Shader<EVAL>(Bling.Shaders.BaseShaderEval<EVAL> txt, Bling.Ops.Operation<IVec<double, Dim<ROW, COL>>, IVec<double, Dim<COL, MID>>, IVec<double, Dim<ROW, MID>>> op) {
        return (a,b) => "mul(" + a + ", " + b + ")";
      }
      public override Expr<IVec<double, Dim<ROW, MID>>> Expand(Expr<IVec<double, Dim<ROW, COL>>> argA, Expr<IVec<double, Dim<COL, MID>>> argB) {
        var RowCol = ToBrandD(argA);
        var ColMid = MatrixBl<COL, MID>.ToBrandD(argB);

        var Es = new DoubleBl[Dim<ROW>.Count,Dim<MID>.Count];
        for (int row = 0; row < Dim<ROW>.Count; row++)
          for (int mid = 0; mid < Dim<MID>.Count; mid++) {
            Es[row, mid] = RowCol.Rows[row].Dot(ColMid.Columns[mid]);
          }
        return new MatrixBl<ROW, MID>(Es).Provider;
      }
    }
    public static MatrixBl<ROW, D1> operator *(BaseMatrixBl<ROW, COL, BRAND> opA, MatrixBl<COL, D1> opB) {
      return opA.Product<D1>(opB); 
    }
    public static MatrixBl<D1,COL> operator *(MatrixBl<D1,ROW> opA, BaseMatrixBl<ROW, COL, BRAND> opB) {
      return opA.Product<COL>(opB);
    }
    public static MatrixBl<ROW, ROW> operator *(BaseMatrixBl<ROW, COL, BRAND> opA, MatrixBl<COL, ROW> opB) {
      return opA.Product<ROW>(opB);
    }
    public MatrixBl<COL, ROW> Transpose {
      get { return TransposeOp.Instance.Make(Provider); }
    }
  }
  public abstract partial class BaseBaseMatrixBl<ROW, COL, BRAND> : BaseMatrixBl<ROW, COL, BRAND>
    where ROW : Dim1D<ROW>
    where COL : Dim1D<COL>
    where BRAND : BaseBaseMatrixBl<ROW, COL, BRAND> {
    public BaseBaseMatrixBl(Expr<IVec<double, Dim<ROW, COL>>> Underlying) : base(Underlying) { }
    public BaseBaseMatrixBl(DoubleBl[] Es) : base(Es) { }
    public BaseBaseMatrixBl(DoubleBl[,] Es) : base(Es) { }
    public BaseBaseMatrixBl() : base() { }
    public BaseBaseMatrixBl(Func<int, int, DoubleBl> F) : base((F)) { }
    public BaseBaseMatrixBl(params MatrixBl<D1, COL>[] Rows) : base((Rows)) { }
    public BaseBaseMatrixBl(params MatrixBl<ROW, D1>[] Columns) : base((Columns)) { }
    public static MatrixBl<COL, COL> operator *(MatrixBl<COL, ROW> opA, BaseBaseMatrixBl<ROW, COL, BRAND> opB) {
      return opA.Product<COL>(opB);
    }
  }
  public partial class MatrixBl<ROW, COL> : BaseBaseMatrixBl<ROW, COL, MatrixBl<ROW, COL>> where ROW : Dim1D<ROW> where COL : Dim1D<COL> {
    public MatrixBl(Expr<IVec<double, Dim<ROW, COL>>> Underlying) : base(Underlying) { }
    public MatrixBl(params DoubleBl[] Es) : base(Es) { }
    public MatrixBl(DoubleBl[,] Es) : base(Es) { }
    public MatrixBl() : base() { }
    public MatrixBl(Func<int, int, DoubleBl> F) : base(F) { }
    public MatrixBl(params MatrixBl<D1, COL>[] Rows) : base(Rows) { }
    public MatrixBl(params MatrixBl<ROW, D1>[] Columns) : base(Columns) { }
    public static implicit operator MatrixBl<ROW, COL>(Expr<IVec<double, Dim<ROW, COL>>> M) { return new MatrixBl<ROW, COL>(M); }
  }
  /// <summary>
  /// Tags an object that can be converted into an affine transformation matrix.
  /// </summary>
  /// <typeparam name="MATRIX">Kind of matrix that conversion can occur to (Matrix2Bl, Matrix3Bl, or Matrix4Bl).</typeparam>
  public interface ICanBeAffineMatrix<MATRIX> where MATRIX : ISquareMatrix<MATRIX> {
    MATRIX Matrix { get; }
  }
  public interface ISquareMatrix<BRAND> : IMatrix<BRAND> where BRAND : ISquareMatrix<BRAND> { }
  public interface ISquareMatrix<DIM, BRAND> : IMatrix<DIM, DIM, BRAND>, ISquareMatrix<BRAND>
    where DIM : Dim1D<DIM>
    where BRAND : ISquareMatrix<DIM, BRAND> { }
  public abstract partial class BaseSquareMatrixBl<DIM, BRAND, PREV> : BaseBaseMatrixBl<DIM, DIM, BRAND>, ISquareMatrix<DIM,BRAND>
    where DIM : Dim1D<DIM>, IPrevDim<DIM,PREV>
    where PREV : Dim1D<PREV>, INextDim<PREV,DIM>
    where BRAND : SquareMatrixBl<DIM, BRAND, PREV> {
    public BaseSquareMatrixBl(Expr<IVec<double, Dim<DIM, DIM>>> Underlying) : base(Underlying) { }
    public BaseSquareMatrixBl(params DoubleBl[] Es) : base(Es) { }
    public BaseSquareMatrixBl(DoubleBl[,] Es) : base(Es) { }
    public BaseSquareMatrixBl(params MatrixBl<D1, DIM>[] rows) : base(rows) { }
    public BaseSquareMatrixBl(params MatrixBl<DIM, D1>[] cols) : base(cols) { }
    public BaseSquareMatrixBl() : base() { }
    public BaseSquareMatrixBl(Func<int, int, DoubleBl> F) : base(F) { }

    static BaseSquareMatrixBl() {
      //Default = Identity;
    }

    private static BRAND Identity0;
    public static BRAND Identity {
      get {
        if (((object)Identity0) == null) {
          Identity0 = ToBrandF((row, col) => row == col ? 1 : 0).VerifyAffine;
        }
        return Identity0;
      }
    }

    protected class ProductOp :
      Ops.Operator<IVec<double, Dim<DIM, DIM>>, IVec<double, Dim<DIM, DIM>>, IVec<double, Dim<DIM, DIM>>, ProductOp> {

      public new static ProductOp Instance = new ProductOp();
      private ProductOp() { }
      public override SmartString Format(SmartString argA, SmartString argB) {
        return ProductOp<DIM>.Instance.Format(argA, argB);
      }
      public override Expr<IVec<double, Dim<DIM, DIM>>> Expand(Expr<IVec<double, Dim<DIM, DIM>>> argA, Expr<IVec<double, Dim<DIM, DIM>>> argB) {
        return ProductOp<DIM>.Instance.Make(argA, argB);
      }
      public override Expr<IVec<double, Dim<DIM, DIM>>> InverseA(Expr<IVec<double, Dim<DIM, DIM>>> argA, Expr<IVec<double, Dim<DIM, DIM>>> argB, Expr<IVec<double, Dim<DIM, DIM>>> result) {
        var M = ToBrand(result);
        var B = ToBrand(argB);
        return M * B.Inverse;
      }
      public override Expr<IVec<double, Dim<DIM, DIM>>> InverseB(Expr<IVec<double, Dim<DIM, DIM>>> argA, Expr<IVec<double, Dim<DIM, DIM>>> argB, Expr<IVec<double, Dim<DIM, DIM>>> result) {
        var M = ToBrand(result);
        var A = ToBrand(argA);
        return A.Inverse * M;
      }
    }
    /// <summary>
    /// Square matrix product, preserves affine transformation properites.
    /// </summary>
    /// <param name="opB"></param>
    /// <returns></returns>
    public virtual BRAND Product(BRAND opB) {
      if (Provider.Equals(Identity.Provider)) return opB;
      if (opB.Provider.Equals(Identity.Provider)) return this;
      var ret = ToBrandD(ProductOp.Instance.Make(this, opB));
      if (opB.IsAffine && this.IsAffine) ret = ret.VerifyAffine;
      return ret;
    }
    public static BRAND operator *(BaseSquareMatrixBl<DIM, BRAND, PREV> opA, BRAND opB) { return opA.Product(opB); }
    public static BRAND operator *(BaseSquareMatrixBl<DIM, BRAND, PREV> opA, ICanBeAffineMatrix<BRAND> opB) { 
      return opA.Product(opB.Matrix.VerifyAffine); 
    }


    public class GenericDetOp : UnaryOperatorS<double,GenericDetOp>, Ops.IMemoOperator<double> {
      public new static readonly GenericDetOp Instance = new GenericDetOp();
      private GenericDetOp() { }
      public override string Name {
        get { return "det"; }
      }
      public override Expr<double> Expand(Expr<IVec<double, Dim<DIM, DIM>>> argA) {
        DoubleBl ret = 0;
        (Dim<DIM>.Count >= 1).Assert();
        var Matrix = ToBrandD(argA);
        // choose the last row because it might be 0, 0, 0, 1
        for (int i = 0; i < Dim<DIM>.Count; i++) {
          ret += Matrix[Dim<DIM>.Count - 1, i] * Matrix.Cofactor(Dim<DIM>.Count - 1, i);
        }
        return ret;
      }
    }
    public class AffineDetOp : UnaryOperatorS<double, AffineDetOp>, Ops.IMemoOperator<double> {
      public new static readonly AffineDetOp Instance = new AffineDetOp();
      private AffineDetOp() { }
      public override string Name {
        get { return "det"; }
      }
      public override Expr<double> Expand(Expr<IVec<double, Dim<DIM, DIM>>> argA) {
        var Matrix = ToBrandD(argA);
        // last row is 0, 0, 0, 1
        return Matrix.Cofactor(Dim<DIM>.Count - 1, Dim<DIM>.Count - 1); 
      }
    }
    /// <summary>
    /// Matrix determanent.
    /// </summary>
    public virtual DoubleBl Det {
      get {
        if (IsAffine) return AffineDetOp.Instance.Make(this);
        return GenericDetOp.Instance.Make(this);
      }
    }
    /// <summary>
    /// Is this matrix an affine transformation matrix?
    /// </summary>
    public bool IsAffine {
      get;
      private set;
    }
    public override BRAND CopyMetaTo(BRAND Other) {
      Other = base.CopyMetaTo(Other);
      Other.IsAffine = this.IsAffine;
      return Other;
    }


    public BoolBl IsReallyAffine {
      get {
        if (IsAffine) return true;
        BoolBl ret = true;
        for (int i = 0; i < Dim<DIM>.Count - 1; i++)
          ret = ret & (M[Dim<DIM>.Count - 1, i] == (i == Dim<DIM>.Count ? 1 : 0));
        return ret;
      }
    }
    public BRAND VerifyAffine {
      get {
        if (this.IsAffine) return (BRAND) this;
        this.IsAffine = true;
        var OldVerifyRhs = this.VerifyRHS;
        this.VerifyRHS = (value) => value.IsAffine && OldVerifyRhs(value);
        this.RHSFail += "; not affine";
        return (BRAND) this;
      }
    }

    public class GenericInverseOp : UnaryOp<GenericInverseOp>, Ops.IMemoOperator<IVec<double,Dim<DIM,DIM>>> {
      public new static readonly GenericInverseOp Instance = new GenericInverseOp();
      private GenericInverseOp() { }
      public override string Name {
        get { return "inverse"; }
      }
      public override Func<Expr<double>, Expr<double>> Uniform {
        get { return null; }
      }
      protected override Expr<IVec<double, Dim<DIM, DIM>>> Expand0(Expr<IVec<double, Dim<DIM, DIM>>> argA) {
        var ArgA = ToBrandD(argA);
        return ToBrandF((row, col) => ArgA.Cofactor(row, col)).Transpose / ArgA.Det;
      }
      public override Expr<IVec<double, Dim<DIM, DIM>>> Inverse(Expr<IVec<double, Dim<DIM, DIM>>> argA, Expr<IVec<double, Dim<DIM, DIM>>> result) {
        return ToBrandD(result).Inverse;
      }
    }
    public class AffineInverseOp : UnaryOp<AffineInverseOp>, Ops.IMemoOperator<IVec<double, Dim<DIM, DIM>>> {
      public new static readonly AffineInverseOp Instance = new AffineInverseOp();
      private AffineInverseOp() { }
      public override string Name {
        get { return "inverse"; }
      }
      public override Func<Expr<double>, Expr<double>> Uniform {
        get { return null; }
      }
      protected override Expr<IVec<double, Dim<DIM, DIM>>> Expand0(Expr<IVec<double, Dim<DIM, DIM>>> argA) {
        var ArgA = ToBrandD(argA).VerifyAffine;
        // inverse(A), b * inverse(A)
        var Am = ArgA.AffineInverse;
        var Bm = (Am * ArgA.Translation) * -1d;
        var Cm = Am.Splice<PREV,PREV,MatrixBl<PREV,PREV>,DIM>(Dim<DIM>.Count - 1, Bm);
        //var XX = new MatrixBl<D1, DIM>((row, col) => col == Dim<DIM>.Count - 1 ? 1 : 0);

        var Dm = Cm.Splice<PREV, DIM, MatrixBl<PREV, DIM>, DIM>(Dim<DIM>.Count - 1, ArgA.Rows[Dim<DIM>.Count-1]);

        return Dm;
      }
      public override Expr<IVec<double, Dim<DIM, DIM>>> Inverse(Expr<IVec<double, Dim<DIM, DIM>>> argA, Expr<IVec<double, Dim<DIM, DIM>>> result) {
        return ToBrandD(result).VerifyAffine.Inverse;
      }
    }
    /// <summary>
    /// Inverse of this matrix.
    /// </summary>
    public virtual BRAND Inverse {
      get {
        if (IsAffine) {
          return ToBrandD(AffineInverseOp.Instance.Make(this)).VerifyAffine;
        }  else return ToBrandD(GenericInverseOp.Instance.Make(this));
      }
    }
    public DoubleBl Minor(int Row, int Column) {
      var matrix = this.DropColumn<DIM, DIM, BRAND, PREV>(Column).DropRow<DIM, PREV, MatrixBl<DIM, PREV>, PREV>(Row);
      return DetP(matrix);
    }
    protected abstract DoubleBl DetP(MatrixBl<PREV, PREV> M);
    /// <summary>
    /// Co-factor of matrix for specified row and column.
    /// </summary>
    public DoubleBl Cofactor(int Row, int Column) {
      return Minor(Row, Column) * (((Row + Column) % 2 == 0).Condition(+1.Bl(),-1));
    }
    /// <summary>
    /// The inverse of the affine section of this affine transformation matrix.
    /// </summary>
    public abstract MatrixBl<PREV, PREV> AffineInverse { get; }

    /// <summary>
    /// Return the northwest submatrix of this affine transformation matrix (the result is of course not affine).
    /// </summary>
    public MatrixBl<PREV, PREV> Affine {
      get {
        IsAffine.Assert();
        var v1 = this.DropColumn<DIM, DIM, BRAND, PREV>(Dim<DIM>.Count - 1);
        var v2 = v1.DropRow<DIM, PREV, MatrixBl<DIM,PREV>, PREV>(Dim<DIM>.Count - 1);
        return v2;
      }
    }
    /// <summary>
    /// Return the translation column of this affine transformation matrix.
    /// </summary>
    public MatrixBl<PREV,D1> Translation {
      get {
        IsAffine.Assert();
        var v1 = this.Columns[Dim<DIM>.Count - 1];
        var v2 = v1.DropRow<DIM, D1, MatrixBl<DIM, D1>, PREV>(Dim<DIM>.Count - 1);
        return v2;
      }
    }

  }
  public abstract partial class SquareMatrixBl<DIM, BRAND, PREV> : BaseSquareMatrixBl<DIM, BRAND, PREV>
    where DIM : Dim1D<DIM>, IPrevDim<DIM, PREV>
    where PREV : Dim1D<PREV>, INextDim<PREV, DIM>
    where BRAND : SquareMatrixBl<DIM, BRAND, PREV> {
    public SquareMatrixBl(Expr<IVec<double, Dim<DIM, DIM>>> Underlying) : base(Underlying) { }
    public SquareMatrixBl(params DoubleBl[] Es) : base(Es) { }
    public SquareMatrixBl(DoubleBl[,] Es) : base(Es) { }
    public SquareMatrixBl(params MatrixBl<D1, DIM>[] rows) : base(rows) { }
    public SquareMatrixBl(params MatrixBl<DIM, D1>[] cols) : base(cols) { }
    public SquareMatrixBl() : base() { }
    public SquareMatrixBl(Func<int, int, DoubleBl> F) : base(F) { }
    /// <summary>
    /// Multiply two square matrices together.
    /// </summary>
    public static BRAND operator *(BRAND opA, SquareMatrixBl<DIM, BRAND, PREV>  opB) {
      return opA.Product((BRAND) opB);
    }
    /// <summary>
    /// Multiply two square matrices together.
    /// </summary>
    public static BRAND operator *(ICanBeAffineMatrix<BRAND> opA, SquareMatrixBl<DIM, BRAND, PREV> opB) {
      return opA.Matrix.VerifyAffine.Product((BRAND)opB);
    }
  }
  public abstract partial class SquareMatrixBl<DIM, BRAND, PREV, POINT> : SquareMatrixBl<DIM, BRAND, PREV>
    where DIM : Dim1D<DIM>, IPrevDim<DIM, PREV>
    where PREV : Dim1D<PREV>, INextDim<PREV, DIM>
    where BRAND : SquareMatrixBl<DIM, BRAND, PREV, POINT>
    where POINT : FixDoubleBl<PREV,POINT> {
    public SquareMatrixBl(Expr<IVec<double, Dim<DIM, DIM>>> Underlying) : base(Underlying) { }
    public SquareMatrixBl(params DoubleBl[] Es) : base(Es) { }
    public SquareMatrixBl(DoubleBl[,] Es) : base(Es) { }
    public SquareMatrixBl(params MatrixBl<D1, DIM>[] rows) : base(rows) { }
    public SquareMatrixBl(params MatrixBl<DIM, D1>[] cols) : base(cols) { }
    public SquareMatrixBl() : base() { }
    public SquareMatrixBl(Func<int, int, DoubleBl> F) : base(F) { }
    /// <summary>
    /// Create a translation affine transform matrix.
    /// </summary>
    /// <param name="By">Translation vector.</param>
    /// <returns>An affine translation matrix.</returns>
    public static BRAND MakeTranslate(POINT By) {
      if (By.Provider.Equals(FixDoubleBl<PREV, POINT>.Zero.Provider)) return Identity;
      var ret = ToBrandF((row, col) => {
        if (row < Dim<PREV>.Count && col == Dim<PREV>.Count ) return By[row];
        else return Identity[row, col];
      });
      return ret.VerifyAffine;
    }
    /// <summary>
    /// Applies a translation affine transform matrix after the current one (translate * this).
    /// </summary>
    /// <param name="By">Translation vector.</param>
    /// <returns>An affine matrix with a translation applied.</returns>
    public BRAND Translate(POINT By) {
      this.IsAffine.Assert();
      return MakeTranslate(By) * this;
    }

    /// <summary>
    /// Create a scale affine transform matrix that scales around a specified center.
    /// </summary>
    /// <param name="Scale">Scale vector matching the dimensions of the desired matrix.</param>
    /// <param name="Center">Center around which to scale (can be zero).</param>
    /// <returns>An affine scale matrix.</returns>
    public static BRAND MakeScale(POINT Scale, POINT Center) {
      if (Scale.Provider.Equals(FixDoubleBl<PREV, POINT>.One.Provider)) return Identity;
      var ret = ToBrandF((row, col) => {
        if (col == row && col < Dim<PREV>.Count) return Scale[col];
        else if (row < Dim<DIM>.Count - 1 && col != row && col == Dim<DIM>.Count - 1) return Center[row] - Scale[row] * Center[row];
        else return Identity[row, col];
      });
      return ret.VerifyAffine;
    }
    /// <summary>
    /// Composes a scale affine transform matrix that scales around a specified center.
    /// </summary>
    /// <param name="Scale">Scale vector matching the dimensions of the desired matrix.</param>
    /// <param name="Center">Center around which to scale (can be zero).</param>
    public BRAND Scale(POINT Scale, POINT Center) {
      this.IsAffine.Assert();
      return MakeScale(Scale, Center) * this;
    }

    /// <summary>
    /// Create a scale affine transform matrix that scales around the origin.
    /// </summary>
    /// <param name="Scale">Scale vector matching the dimensions of the desired matrix.</param>
    /// <returns>An affine scale matrix.</returns>
    public static BRAND MakeScale(POINT Scale0) { return MakeScale(Scale0, FixDoubleBl<PREV, POINT>.Repeat(0d)); }
    /// <summary>
    /// Applies a scale affine transform matrix that scales around the origin.
    /// </summary>
    /// <param name="Scale">Scale vector matching the dimensions of the desired matrix.</param>
    public BRAND Scale(POINT Scale) {
      this.IsAffine.Assert();
      return MakeScale(Scale) * this;
    }


    /// <summary>
    /// Transform point according to this affine transformation matrix.
    /// </summary>
    public POINT Transform(POINT P) {
      this.IsAffine.Assert();
      MatrixBl<PREV,D1> M0 = P;
      var M1 = M0.Splice<PREV, D1, MatrixBl<PREV, D1>, DIM>(Dim<PREV>.Count, 1d.Bl());
      var M2 = this * M1;
      var last = M2[Dim<PREV>.Count, 0];
      var M3 = M2.DropRow<DIM, D1, MatrixBl<DIM, D1>, PREV>(Dim<PREV>.Count);
      M3 = M3 / last;
      return FixDoubleBl<PREV,POINT>.Convert(M3); 
    }
  }
  public partial class Matrix2Bl : SquareMatrixBl<D2, Matrix2Bl, D1, Double1Bl> {
    public Matrix2Bl(Expr<IVec<double, Dim<D2, D2>>> Underlying) : base(Underlying) { }
    public Matrix2Bl(params DoubleBl[] Es) : base(Es) { }
    public Matrix2Bl(DoubleBl[,] Es) : base(Es) { }
    public Matrix2Bl(MatrixBl<D1, D2> rowA, MatrixBl<D1, D2> rowB) : base(rowA, rowB) { }
    public Matrix2Bl(MatrixBl<D2, D1> colA, MatrixBl<D2, D1> colB) : base(colA, colB) { }
    public Matrix2Bl() : base() { }
    public Matrix2Bl(Func<int, int, DoubleBl> F) : base(F) { }
    public static implicit operator MatrixBl<D2, D2>(Matrix2Bl M) { return new MatrixBl<D2, D2>(M.Provider); }
    public static implicit operator Matrix2Bl(MatrixBl<D2, D2> M) { return new Matrix2Bl(M.Provider); }
    public static implicit operator Matrix2Bl(Expr<IVec<double,Dim<D2,D2>>> M) { return new Matrix2Bl(M); }
    public static Matrix2Bl ByColumns(Double2Bl rowA, Double2Bl rowB) {
      return new Matrix2Bl((MatrixBl<D2, D1>)rowA, (MatrixBl<D2, D1>)rowB);
    }
    public static Matrix2Bl ByRows(Double2Bl colA, Double2Bl colB) {
      return new Matrix2Bl((MatrixBl<D1, D2>)colA, (MatrixBl<D1, D2>)colB);
    }
    //public override DoubleBl Det { get { return this[0, 0] * this[1, 1] - this[0, 1] * this[1, 0]; } }
    protected override DoubleBl DetP(MatrixBl<D1, D1> M) { return M.Det(); }
    public override MatrixBl<D1, D1> AffineInverse { get { return Affine; } }
  }

  public partial class Matrix3Bl : SquareMatrixBl<D3, Matrix3Bl, D2, Double2Bl> {
    public Matrix3Bl(Expr<IVec<double, Dim<D3, D3>>> Underlying) : base(Underlying) { }
    public Matrix3Bl(params DoubleBl[] Es) : base(Es) { }
    public Matrix3Bl(DoubleBl[,] Es) : base(Es) { }
    public Matrix3Bl(MatrixBl<D1, D3> rowA, MatrixBl<D1, D3> rowB, MatrixBl<D1, D3> rowC) : base(rowA, rowB, rowC) { }
    public Matrix3Bl(MatrixBl<D3, D1> colA, MatrixBl<D3, D1> colB, MatrixBl<D3, D1> colC) : base(colA, colB, colC) { }
    public Matrix3Bl() : base() { }
    public Matrix3Bl(Func<int, int, DoubleBl> F) : base(F) { }
    public static implicit operator MatrixBl<D3, D3>(Matrix3Bl M) { return new MatrixBl<D3, D3>(M.Provider); }
    public static implicit operator Matrix3Bl(MatrixBl<D3, D3> M) { return new Matrix3Bl(M.Provider); }
    public static implicit operator Matrix3Bl(Expr<IVec<double, Dim<D3, D3>>> M) { return new Matrix3Bl(M); }
    public static Matrix3Bl ByColumns(Double3Bl colA, Double3Bl colB, Double3Bl colC) {
      return new Matrix3Bl((MatrixBl<D3, D1>)colA, (MatrixBl<D3, D1>)colB, (MatrixBl<D3, D1>)colC);
    }
    public static Matrix3Bl ByRows(Double3Bl rowA, Double3Bl rowB, Double3Bl rowC) {
      return new Matrix3Bl((MatrixBl<D1, D3>)rowA, (MatrixBl<D1, D3>)rowB, (MatrixBl<D1, D3>)rowC);
    }
    protected override DoubleBl DetP(MatrixBl<D2, D2> M) { return M.Det(); }
    public override MatrixBl<D2, D2> AffineInverse { get { return Affine.Bl().Inverse; } }
  }
  public partial class Matrix4Bl : SquareMatrixBl<D4, Matrix4Bl, D3, Double3Bl> {
    public Matrix4Bl(Expr<IVec<double, Dim<D4, D4>>> Underlying) : base(Underlying) { }
    public Matrix4Bl(params DoubleBl[] Es) : base(Es) { }
    public Matrix4Bl(DoubleBl[,] Es) : base(Es) { }
    public Matrix4Bl(MatrixBl<D1, D4> rowA, MatrixBl<D1, D4> rowB, MatrixBl<D1, D4> rowC, MatrixBl<D1, D4> rowD) : base(rowA, rowB, rowC, rowD) { }
    public Matrix4Bl(MatrixBl<D4, D1> colA, MatrixBl<D4, D1> colB, MatrixBl<D4, D1> colC, MatrixBl<D4, D1> colD) : base(colA, colB, colC, colD) { }
    public Matrix4Bl() : base() { }
    public Matrix4Bl(Func<int, int, DoubleBl> F) : base(F) { }
    public static implicit operator MatrixBl<D4, D4>(Matrix4Bl M) { return new MatrixBl<D4, D4>(M.Provider); }
    public static implicit operator Matrix4Bl(MatrixBl<D4, D4> M) { return new Matrix4Bl(M.Provider); }
    public static implicit operator Matrix4Bl(Expr<IVec<double, Dim<D4, D4>>> M) { return new Matrix4Bl(M); }
    public static Matrix4Bl ByColumns(Double4Bl colA, Double4Bl colB, Double4Bl colC, Double4Bl colD) {
      return new Matrix4Bl((MatrixBl<D4, D1>)colA, (MatrixBl<D4, D1>)colB, (MatrixBl<D4, D1>)colC, (MatrixBl<D4, D1>)colD);
    }
    public static Matrix4Bl ByRows(Double4Bl rowA, Double4Bl rowB, Double4Bl rowC, Double4Bl rowD) {
      return new Matrix4Bl((MatrixBl<D1, D4>)rowA, (MatrixBl<D1, D4>)rowB, (MatrixBl<D1, D4>)rowC, (MatrixBl<D1, D4>)rowD);
    }
    protected override DoubleBl DetP(MatrixBl<D3, D3> M) { return M.Det(); }
    //public Double3Bl Transform(Double3Bl p) { return (Double3Bl)Transform((Double4Bl)p); }
    public override MatrixBl<D3, D3> AffineInverse { get { return Affine.Bl().Inverse; } }
  }
  public static class MatrixExtensions {

    public static FlexDoubleBl<D> Bl<D>(this MatrixBl<D1, D> M) where D : Dim1D<D> { return VecMatrix<D>.FromMatrixA.Instance.Make(M); }
    public static FlexDoubleBl<D> Bl<D>(this MatrixBl<D, D1> M) where D : Dim1D<D> { return VecMatrix<D>.FromMatrixA.Instance.Make(M.Transpose); }

    public static DoubleBl Dot<D>(this MatrixBl<D1, D> mA, MatrixBl<D, D1> mB)
      where D : Dim1D<D> { return mA.Bl().Dot(mB.Bl()); }
    public static DoubleBl Dot<D>(this MatrixBl<D, D1> mB, MatrixBl<D1, D> mA)
      where D : Dim1D<D> { return mA.Bl().Dot(mB.Bl()); }

    public static MatrixBl<ROW_NEXT, COL> Splice<ROW, COL, BRAND, ROW_NEXT>(this BaseBaseMatrixBl<ROW, COL, BRAND> Matrix, int Index, MatrixBl<D1, COL> Row)
      where ROW : Dim1D<ROW>, INextDim<ROW, ROW_NEXT>
      where COL : Dim1D<COL>
      where ROW_NEXT : Dim1D<ROW_NEXT>, IPrevDim<ROW_NEXT, ROW>
      where BRAND : BaseBaseMatrixBl<ROW, COL, BRAND> {
      if (typeof(COL) == typeof(D1)) {
        var IndicesA = new int[Index];
        for (int i = 0; i < Index; i++) IndicesA[i] = i;
        var IndicesB = new int[Dim<ROW>.Count - Index];
        for (int i = Index; i < Dim<ROW>.Count;  i++) IndicesB[i] = i;
        if (IndicesB.Length == 0) {
          // this is just a combine
          var U = MatrixBl<ROW_NEXT, COL>.CombineOp<Dim<ROW, COL>, Dim<D1, COL>>.Instance.Make(Matrix.UnderlyingD, Row.UnderlyingD);
          return new MatrixBl<ROW_NEXT, COL>(U);
        } else if (IndicesA.Length == 0) {
          var U = MatrixBl<ROW_NEXT, COL>.CombineOp<Dim<D1, COL>, Dim<ROW, COL>>.Instance.Make(Row.UnderlyingD, Matrix.UnderlyingD);
          return new MatrixBl<ROW_NEXT, COL>(U);
        } 
         




        true.Assert();
      }

      return new MatrixBl<ROW_NEXT, COL>((row, col) => 
        (row < Index) ? Matrix[row, col] :
        (row == Index) ? Row[0, col] : Matrix[row - 1, col]);
    }
    public static MatrixBl<D3, COL> Splice<COL, BRAND>(this BaseBaseMatrixBl<D2, COL, BRAND> Matrix, int Index, MatrixBl<D1, COL> Row)
      where COL : Dim1D<COL>
      where BRAND : BaseBaseMatrixBl<D2, COL, BRAND> {
      return Splice<D2, COL, BRAND, D3>(Matrix, Index, Row);
    }
    public static MatrixBl<D4, COL> Splice<COL, BRAND>(this BaseBaseMatrixBl<D3, COL, BRAND> Matrix, int Index, MatrixBl<D1, COL> Row)
      where COL : Dim1D<COL>
      where BRAND : BaseBaseMatrixBl<D3, COL, BRAND> {
      return Splice<D3, COL, BRAND, D4>(Matrix, Index, Row);
    }
    public static MatrixBl<ROW, COL_NEXT> Splice<ROW, COL, BRAND, COL_NEXT>(this BaseBaseMatrixBl<ROW, COL, BRAND> Matrix, int Index, MatrixBl<ROW, D1> Row)
      where ROW : Dim1D<ROW>
      where COL : Dim1D<COL>, INextDim<COL, COL_NEXT>
      where COL_NEXT : Dim1D<COL_NEXT>, IPrevDim<COL_NEXT, COL>
      where BRAND : BaseBaseMatrixBl<ROW, COL, BRAND> {
      if (typeof(ROW) == typeof(D1)) {
        true.Assert();
      }

      return new MatrixBl<ROW, COL_NEXT>((row, col) => 
        (col < Index) ? Matrix[row, col] :
        (col == Index) ? Row[row, 0] : Matrix[row, col - 1]);
    }
    public static MatrixBl<ROW, D3> Splice<ROW, BRAND>(this BaseBaseMatrixBl<ROW, D2, BRAND> Matrix, int Index, MatrixBl<ROW, D1> Row)
      where ROW : Dim1D<ROW>
      where BRAND : BaseBaseMatrixBl<ROW, D2, BRAND> {
      return Splice<ROW, D2, BRAND, D3>(Matrix, Index, Row);
    }
    public static MatrixBl<ROW, D4> Splice<ROW, BRAND>(this BaseBaseMatrixBl<ROW, D3, BRAND> Matrix, int Index, MatrixBl<ROW, D1> Row)
      where ROW : Dim1D<ROW>
      where BRAND : BaseBaseMatrixBl<ROW, D3, BRAND> {
      return Splice<ROW, D3, BRAND, D4>(Matrix, Index, Row);
    }


    public static MatrixBl<ROW_PREV, COL> DropRow<ROW, COL, BRAND, ROW_PREV>(this BaseBaseMatrixBl<ROW, COL, BRAND> Matrix, int Index)
      where ROW : Dim1D<ROW>, IPrevDim<ROW, ROW_PREV>
      where COL : Dim1D<COL>
      where ROW_PREV : Dim1D<ROW_PREV>, INextDim<ROW_PREV, ROW>
      where BRAND : BaseBaseMatrixBl<ROW, COL, BRAND> {
      if (typeof(COL) == typeof(D1)) {
        // can we try to swizzle?
        var Indices = new int[Dim<ROW_PREV>.Count];
        for (int i = 0; i < Index; i++) Indices[i] = i;
        for (int i = Index + 1; i < Dim<ROW>.Count; i++) Indices[i - 1] = i;
        var Next = Matrix.Swizzle<Dim<ROW_PREV, COL>>(Indices);

        return new MatrixBl<ROW_PREV, COL>(Next.Provider);
      }
      return new MatrixBl<ROW_PREV, COL>((row, col) => (row < Index) ? Matrix[row, col] : Matrix[row + 1, col]);
    }
    public static MatrixBl<D1, COL> DropRow<COL, BRAND>(this BaseBaseMatrixBl<D2, COL, BRAND> Matrix, int Index)
      where COL : Dim1D<COL>
      where BRAND : BaseBaseMatrixBl<D2, COL, BRAND> {
      return DropRow<D2, COL, BRAND, D1>(Matrix, Index);
    }
    public static MatrixBl<D2, COL> DropRow<COL, BRAND>(this BaseBaseMatrixBl<D3, COL, BRAND> Matrix, int Index)
      where COL : Dim1D<COL>
      where BRAND : BaseBaseMatrixBl<D3, COL, BRAND> {
      return DropRow<D3, COL, BRAND, D2>(Matrix, Index);
    }
    public static MatrixBl<D3, COL> DropRow<COL, BRAND>(this BaseBaseMatrixBl<D4, COL, BRAND> Matrix, int Index)
      where COL : Dim1D<COL>
      where BRAND : BaseBaseMatrixBl<D4, COL, BRAND> {
      return DropRow<D4, COL, BRAND, D3>(Matrix, Index);
    }

    public static MatrixBl<ROW, COL_PREV> DropColumn<ROW, COL, BRAND, COL_PREV>(this BaseBaseMatrixBl<ROW, COL, BRAND> Matrix, int Index)
      where ROW : Dim1D<ROW>
      where COL : Dim1D<COL>, IPrevDim<COL, COL_PREV>
      where COL_PREV : Dim1D<COL_PREV>, INextDim<COL_PREV, COL>
      where BRAND : BaseBaseMatrixBl<ROW, COL, BRAND> {
      if (typeof(ROW) == typeof(D1)) {

        true.Assert();
      }

      return new MatrixBl<ROW, COL_PREV>((row, col) => (col < Index) ? Matrix[row, col] : Matrix[row, col + 1]);
    }
    public static MatrixBl<ROW, D1> DropColumn<ROW, BRAND>(this BaseBaseMatrixBl<ROW, D2, BRAND> Matrix, int Index)
      where ROW : Dim1D<ROW>
      where BRAND : BaseBaseMatrixBl<ROW, D2, BRAND> {
      return DropColumn<ROW, D2, BRAND, D1>(Matrix, Index);
    }
    public static MatrixBl<ROW, D2> DropColumn<ROW, BRAND>(this BaseBaseMatrixBl<ROW, D3, BRAND> Matrix, int Index)
      where ROW : Dim1D<ROW>
      where BRAND : BaseBaseMatrixBl<ROW, D3, BRAND> {
      return DropColumn<ROW, D3, BRAND, D2>(Matrix, Index);
    }
    public static MatrixBl<ROW, D3> DropColumn<ROW, BRAND>(this BaseBaseMatrixBl<ROW, D4, BRAND> Matrix, int Index)
      where ROW : Dim1D<ROW>
      where BRAND : BaseBaseMatrixBl<ROW, D4, BRAND> {
      return DropColumn<ROW, D4, BRAND, D3>(Matrix, Index);
    }
    public static Matrix3Bl Quardrant(this Matrix4Bl matrix, bool isTop, bool isLeft) {
      return matrix.DropRow(isTop ? 0 : 3).DropColumn(isLeft ? 0 : 3);
    }
    public static Matrix2Bl Quardrant(this Matrix3Bl matrix, bool isTop, bool isLeft) {
      return matrix.DropRow(isTop ? 0 : 2).DropColumn(isLeft ? 0 : 2);
    }
    public static DoubleBl Det(this MatrixBl<D1, D1> M) { return M; }
    public static DoubleBl Det(this MatrixBl<D2, D2> M) { return M.Bl().Det; }
    public static DoubleBl Det(this MatrixBl<D3, D3> M) { return M.Bl().Det; }
    public static DoubleBl Det(this MatrixBl<D4, D4> M) { return M.Bl().Det; }


    public static DoubleBl DetOld(this Matrix4Bl A) {
      return A[0, 0] * A[1, 1] - A[1, 0] * A[0, 1];
    }

  }
}
